// Robustness checks for simulation of emissions pricing
// Author Levi Marks levi.marks.1>at>gmail.com

// Preliminaries
global DATA "[Insert Path Here]"
global iterations 1000
global step_size .05
set more off
set matsize 10000
graph drop _all

local iterations 670
local step_size $step_size
local ch4_cont 0.79

// Preliminary regression macros
local dependent			rate
local independent  		p_spot 
local controls		 	wells oil completions
local FE 				i.region_year i.facility_id 
local conditions		if trim_5 != 1


// 0 // Main specification

local name main

use "$DATA/ghgrp_di_intermediate", clear

scalar dimension = 2
fp <p_spot>, replace dimension(`=dimension') powers(-2 -2): reg `dependent' <p_spot> `controls' `FE'  `conditions', vce(cluster facility)
keep `conditions'

do "$DATA/3_sim_setup.do"

forval i = 1/`iterations' {
	gen price_`i' = price_00 + `ch4_cont'*`step_size'*`i'
	
	// First derivative
	if dimension == 1 & fp_1 != 0 {
		gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' )
	}
	if dimension == 2 & fp_1 != 0 & fp_2 != 0 {
		if (fp_1 != fp_2) {
			gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*fp_2*price_`i'^`=fp_2-1' )
		}
		if (fp_1 == fp_2) {
			gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*( fp_2*(price_`i'^`=fp_2-1')*log(price_`i') + (price_`i'^`=fp_2')*(1/price_`i') ) )
		}
	}
	
	qui replace rate_`i' = lower_bound if rate_`i' < lower_bound
	qui replace rate_`i' = rate_00 if rate_00 < lower_bound 
	gen emissions_`i' = rate_`i'*gas_0
	gen reduction_`i' = emissions_`=(`i'-1)' - emissions_`i'
	gen value_`i' = reduction_`i'*price_00 
	gen cost_`i' = ((price_`i' + price_`=(`i'-1)')/2)*(emissions_`=(`i'-1)' - emissions_`i')
	egen sum_cost_`i' = sum(cost_`i')
	egen sum_reduction_`i' = sum(reduction_`i')
	egen sum_value_`i' = sum(value_`i')
	gen ag_cost_`i' = ag_cost_`=`i'-1' + sum_cost_`i'
	gen ag_reduction_`i' = ag_reduction_`=`i'-1' + sum_reduction_`i'
	gen ag_value_`i' = ag_value_`=`i'-1' + sum_value_`i'
	gen mc_`i' = sum_cost_`i'/sum_reduction_`i'
	drop emissions_`=(`i'-1)' reduction_`i' value_`i' cost_`i' sum_* rate_`=(`i'-1)' price_`=(`i'-1)' //Dropping unnecessary variables to increase speed
}

do "$DATA/3_clean_sim_results.do"

preserve
keep mc_mcf mc_tco2e reduction_percent ag_reduction_tco2e
compress*
save "$DATA/Temp/sim_robust_`name'_graph.dta", replace
restore

//do "$DATA/3_sim_graph_macc.do" // Optional graph

gen check = "`name'"
keep if tax_tco2 == 5 | tax_tco2 == 20 | tax_tco2 == 44.4
keep check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
order check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
compress *
save "$DATA/Temp/sim_robust_`name'_table.dta", replace

*/






// 1 // First-order FP

local name FP1

use "$DATA/ghgrp_di_intermediate", clear

scalar dimension = 1
fp <p_spot>, replace dimension(`=dimension') powers(): reg `dependent' <p_spot> `controls' `FE'  `conditions', vce(cluster facility)
keep `conditions'

do "$DATA/3_sim_setup.do"

forval i = 1/`iterations' {
	gen price_`i' = price_00 + `ch4_cont'*`step_size'*`i'
	
	// First derivative
	if dimension == 1 & fp_1 != 0 {
		gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' )
	}
	if dimension == 2 & fp_1 != 0 & fp_2 != 0 {
		if (fp_1 != fp_2) {
			gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*fp_2*price_`i'^`=fp_2-1' )
		}
		if (fp_1 == fp_2) {
			gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*( fp_2*(price_`i'^`=fp_2-1')*log(price_`i') + (price_`i'^`=fp_2')*(1/price_`i') ) )
		}
	}
	qui replace rate_`i' = lower_bound if rate_`i' < lower_bound
	qui replace rate_`i' = rate_00 if rate_00 < lower_bound 
	gen emissions_`i' = rate_`i'*gas_0
	gen reduction_`i' = emissions_`=(`i'-1)' - emissions_`i'
	gen value_`i' = reduction_`i'*price_00 
	gen cost_`i' = ((price_`i' + price_`=(`i'-1)')/2)*(emissions_`=(`i'-1)' - emissions_`i')
	egen sum_cost_`i' = sum(cost_`i')
	egen sum_reduction_`i' = sum(reduction_`i')
	egen sum_value_`i' = sum(value_`i')
	gen ag_cost_`i' = ag_cost_`=`i'-1' + sum_cost_`i'
	gen ag_reduction_`i' = ag_reduction_`=`i'-1' + sum_reduction_`i'
	gen ag_value_`i' = ag_value_`=`i'-1' + sum_value_`i'
	gen mc_`i' = sum_cost_`i'/sum_reduction_`i'
	drop emissions_`=(`i'-1)' reduction_`i' value_`i' cost_`i' sum_* rate_`=(`i'-1)' price_`=(`i'-1)' //Dropping unnecessary variables to increase speed
}

do "$DATA/3_clean_sim_results.do"

preserve
keep mc_mcf mc_tco2e reduction_percent ag_reduction_tco2e
compress*
save "$DATA/Temp/sim_robust_`name'_graph.dta", replace
restore

//do "$DATA/3_sim_graph_macc.do" // Optional graph

gen check = "`name'"
keep if tax_tco2 == 5 | tax_tco2 == 20 | tax_tco2 == 44.4
keep check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
order check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
compress *
save "$DATA/Temp/sim_robust_`name'_table.dta", replace

*/







// 2 // Weighted by average gas production

use "$DATA/ghgrp_di_intermediate", clear

local name weighted

scalar dimension = 2
fp <p_spot>, replace dimension(`=dimension') powers(): reg `dependent' <p_spot> `controls' `FE'  `conditions' [aweight=mean_gas], vce(cluster facility)
keep `conditions'

do "$DATA/3_sim_setup.do"

forval i = 1/`iterations' {
	gen price_`i' = price_00 + `ch4_cont'*`step_size'*`i'
	
	// First derivative
	if dimension == 1 & fp_1 != 0 {
		gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' )
	}
	if dimension == 2 & fp_1 != 0 & fp_2 != 0 {
		if (fp_1 != fp_2) {
			gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*fp_2*price_`i'^`=fp_2-1' )
		}
		if (fp_1 == fp_2) {
			gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*( fp_2*(price_`i'^`=fp_2-1')*log(price_`i') + (price_`i'^`=fp_2')*(1/price_`i') ) )
		}
	}
	qui replace rate_`i' = lower_bound if rate_`i' < lower_bound
	qui replace rate_`i' = rate_00 if rate_00 < lower_bound 
	gen emissions_`i' = rate_`i'*gas_0
	gen reduction_`i' = emissions_`=(`i'-1)' - emissions_`i'
	gen value_`i' = reduction_`i'*price_00 
	gen cost_`i' = ((price_`i' + price_`=(`i'-1)')/2)*(emissions_`=(`i'-1)' - emissions_`i')
	egen sum_cost_`i' = sum(cost_`i')
	egen sum_reduction_`i' = sum(reduction_`i')
	egen sum_value_`i' = sum(value_`i')
	gen ag_cost_`i' = ag_cost_`=`i'-1' + sum_cost_`i'
	gen ag_reduction_`i' = ag_reduction_`=`i'-1' + sum_reduction_`i'
	gen ag_value_`i' = ag_value_`=`i'-1' + sum_value_`i'
	gen mc_`i' = sum_cost_`i'/sum_reduction_`i'
	drop emissions_`=(`i'-1)' reduction_`i' value_`i' cost_`i' sum_* rate_`=(`i'-1)' price_`=(`i'-1)' //Dropping unnecessary variables to increase speed
}

do "$DATA/3_clean_sim_results.do"

preserve
keep mc_mcf mc_tco2e reduction_percent ag_reduction_tco2e
compress*
save "$DATA/Temp/sim_robust_`name'_graph.dta", replace
restore

// do "$DATA/3_sim_graph_macc.do" // Optional graph

gen check = "`name'"
keep if tax_tco2 == 5 | tax_tco2 == 20 | tax_tco2 == 44.4
keep check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
order check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
compress *
save "$DATA/Temp/sim_robust_`name'_table.dta", replace

*/







// 3 // Trimming emission rates at 1%

local name trim1

use "$DATA/ghgrp_di_intermediate", clear

local conditions if trim_1 != 1

scalar dimension = 2
fp <p_spot>, replace dimension(`=dimension') powers(): reg `dependent' <p_spot> `controls' `FE'  `conditions', vce(cluster facility)
keep `conditions'

do "$DATA/3_sim_setup.do"

forval i = 1/`iterations' {
	gen price_`i' = price_00 + `ch4_cont'*`step_size'*`i'
	
	// First derivative
	if dimension == 1 & fp_1 != 0 {
		gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' )
	}
	if dimension == 2 & fp_1 != 0 & fp_2 != 0 {
		if (fp_1 != fp_2) {
			gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*fp_2*price_`i'^`=fp_2-1' )
		}
		if (fp_1 == fp_2) {
			gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*( fp_2*(price_`i'^`=fp_2-1')*log(price_`i') + (price_`i'^`=fp_2')*(1/price_`i') ) )
		}
	}
	qui replace rate_`i' = lower_bound if rate_`i' < lower_bound
	qui replace rate_`i' = rate_00 if rate_00 < lower_bound 
	gen emissions_`i' = rate_`i'*gas_0
	gen reduction_`i' = emissions_`=(`i'-1)' - emissions_`i'
	gen value_`i' = reduction_`i'*price_00 
	gen cost_`i' = ((price_`i' + price_`=(`i'-1)')/2)*(emissions_`=(`i'-1)' - emissions_`i')
	egen sum_cost_`i' = sum(cost_`i')
	egen sum_reduction_`i' = sum(reduction_`i')
	egen sum_value_`i' = sum(value_`i')
	gen ag_cost_`i' = ag_cost_`=`i'-1' + sum_cost_`i'
	gen ag_reduction_`i' = ag_reduction_`=`i'-1' + sum_reduction_`i'
	gen ag_value_`i' = ag_value_`=`i'-1' + sum_value_`i'
	gen mc_`i' = sum_cost_`i'/sum_reduction_`i'
	drop emissions_`=(`i'-1)' reduction_`i' value_`i' cost_`i' sum_* rate_`=(`i'-1)' price_`=(`i'-1)' //Dropping unnecessary variables to increase speed
}

do "$DATA/3_clean_sim_results.do"

preserve
keep mc_mcf mc_tco2e reduction_percent ag_reduction_tco2e
compress*
save "$DATA/Temp/sim_robust_`name'_graph.dta", replace
restore

// do "$DATA/3_sim_graph_macc.do" // Optional graph

gen check = "`name'"
keep if tax_tco2 == 5 | tax_tco2 == 20 | tax_tco2 == 44.4
keep check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
order check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
compress *
save "$DATA/Temp/sim_robust_`name'_table.dta", replace

local conditions if trim_5 != 1 // Reset trim

*/







// 4 // State regs

local name state_regs

use "$DATA/ghgrp_di_intermediate", clear

drop if state_regs == 1
preserve
bysort facility_id: egen min_rate = min(rate)
bysort facility_id: egen max_rate = max(rate)
collapse min_rate max_rate, by(facility_id)
sort min_rate
gen trim_5 = 1 if _n <= round(_N/20)
sort max_rate
replace trim_5 = 1 if _n > _N - round(_N/20)
drop min max
sort facility_id
save "$DATA/Temp/trim_facilities", replace
restore
merge m:1 facility_id using "$DATA/Temp/trim_facilities"
drop _merge

scalar dimension = 2
fp <p_spot>, replace dimension(`=dimension') powers(): reg `dependent' <p_spot> `controls' `FE'  `conditions', vce(cluster facility)
keep `conditions'

do "$DATA/3_sim_setup.do"

forval i = 1/`iterations' {
	gen price_`i' = price_00 + `ch4_cont'*`step_size'*`i'
	
	// First derivative
	if dimension == 1 & fp_1 != 0 {
		gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' )
	}
	if dimension == 2 & fp_1 != 0 & fp_2 != 0 {
		if (fp_1 != fp_2) {
			gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*fp_2*price_`i'^`=fp_2-1' )
		}
		if (fp_1 == fp_2) {
			gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*( fp_2*(price_`i'^`=fp_2-1')*log(price_`i') + (price_`i'^`=fp_2')*(1/price_`i') ) )
		}
	}
	qui replace rate_`i' = lower_bound if rate_`i' < lower_bound
	qui replace rate_`i' = rate_00 if rate_00 < lower_bound 
	gen emissions_`i' = rate_`i'*gas_0
	gen reduction_`i' = emissions_`=(`i'-1)' - emissions_`i'
	gen value_`i' = reduction_`i'*price_00 
	gen cost_`i' = ((price_`i' + price_`=(`i'-1)')/2)*(emissions_`=(`i'-1)' - emissions_`i')
	egen sum_cost_`i' = sum(cost_`i')
	egen sum_reduction_`i' = sum(reduction_`i')
	egen sum_value_`i' = sum(value_`i')
	gen ag_cost_`i' = ag_cost_`=`i'-1' + sum_cost_`i'
	gen ag_reduction_`i' = ag_reduction_`=`i'-1' + sum_reduction_`i'
	gen ag_value_`i' = ag_value_`=`i'-1' + sum_value_`i'
	gen mc_`i' = sum_cost_`i'/sum_reduction_`i'
	drop emissions_`=(`i'-1)' reduction_`i' value_`i' cost_`i' sum_* rate_`=(`i'-1)' price_`=(`i'-1)' //Dropping unnecessary variables to increase speed
}

do "$DATA/3_clean_sim_results.do"

preserve
keep mc_mcf mc_tco2e reduction_percent ag_reduction_tco2e
compress*
save "$DATA/Temp/sim_robust_`name'_graph.dta", replace
restore

// do "$DATA/3_sim_graph_macc.do" // Optional graph

gen check = "`name'"
keep if tax_tco2 == 5 | tax_tco2 == 20 | tax_tco2 == 44.4
keep check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
order check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
compress *
save "$DATA/Temp/sim_robust_`name'_table.dta", replace

*/






// 5 // Controlling for completions in previous year

local name comp_lastyr

use "$DATA/ghgrp_di_intermediate", clear

drop if missing(ghgrp_comp_lastyr)
preserve
bysort facility_id: egen min_rate = min(rate)
bysort facility_id: egen max_rate = max(rate)
collapse min_rate max_rate, by(facility_id)
sort min_rate
gen trim_5 = 1 if _n <= round(_N/20)
sort max_rate
replace trim_5 = 1 if _n > _N - round(_N/20)
drop min max
sort facility_id
save "$DATA/Temp/trim_facilities", replace
restore
merge m:1 facility_id using "$DATA/Temp/trim_facilities"
drop _merge

local controls wells oil completions ghgrp_comp_lastyr

scalar dimension = 2
fp <p_spot>, replace dimension(`=dimension') powers(): reg `dependent' <p_spot> `controls' `FE'  `conditions', vce(cluster facility)
keep `conditions'

do "$DATA/3_sim_setup.do"

forval i = 1/`iterations' {
	gen price_`i' = price_00 + `ch4_cont'*`step_size'*`i'
	
	// First derivative
	if dimension == 1 & fp_1 != 0 {
		gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' )
	}
	if dimension == 2 & fp_1 != 0 & fp_2 != 0 {
		if (fp_1 != fp_2) {
			gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*fp_2*price_`i'^`=fp_2-1' )
		}
		if (fp_1 == fp_2) {
			gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*( fp_2*(price_`i'^`=fp_2-1')*log(price_`i') + (price_`i'^`=fp_2')*(1/price_`i') ) )
		}
	}
	qui replace rate_`i' = lower_bound if rate_`i' < lower_bound
	qui replace rate_`i' = rate_00 if rate_00 < lower_bound 
	gen emissions_`i' = rate_`i'*gas_0
	gen reduction_`i' = emissions_`=(`i'-1)' - emissions_`i'
	gen value_`i' = reduction_`i'*price_00 
	gen cost_`i' = ((price_`i' + price_`=(`i'-1)')/2)*(emissions_`=(`i'-1)' - emissions_`i')
	egen sum_cost_`i' = sum(cost_`i')
	egen sum_reduction_`i' = sum(reduction_`i')
	egen sum_value_`i' = sum(value_`i')
	gen ag_cost_`i' = ag_cost_`=`i'-1' + sum_cost_`i'
	gen ag_reduction_`i' = ag_reduction_`=`i'-1' + sum_reduction_`i'
	gen ag_value_`i' = ag_value_`=`i'-1' + sum_value_`i'
	gen mc_`i' = sum_cost_`i'/sum_reduction_`i'
	drop emissions_`=(`i'-1)' reduction_`i' value_`i' cost_`i' sum_* rate_`=(`i'-1)' price_`=(`i'-1)' //Dropping unnecessary variables to increase speed
}

do "$DATA/3_clean_sim_results.do"

preserve
keep mc_mcf mc_tco2e reduction_percent ag_reduction_tco2e
compress*
save "$DATA/Temp/sim_robust_`name'_graph.dta", replace
restore

// do "$DATA/3_sim_graph_macc.do" // Optional graph

gen check = "`name'"
keep if tax_tco2 == 5 | tax_tco2 == 20 | tax_tco2 == 44.4
keep check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
order check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
compress *
save "$DATA/Temp/sim_robust_`name'_table.dta", replace

local controls wells oil completions // Reset

*/







// 6 // Year FE

local name year_FE

use "$DATA/ghgrp_di_intermediate", clear

local FE i.year i.facility_id

scalar dimension = 1
fp <p_spot>, replace dimension(`=dimension') powers(): reg `dependent' <p_spot> `controls' `FE'  `conditions', vce(cluster facility)
keep `conditions'

do "$DATA/3_sim_setup.do"

forval i = 1/`iterations' {
	gen price_`i' = price_00 + `ch4_cont'*`step_size'*`i'
	
	// First derivative
	if dimension == 1 & fp_1 != 0 {
		gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' )
	}
	if dimension == 2 & fp_1 != 0 & fp_2 != 0 {
		if (fp_1 != fp_2) {
			gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*fp_2*price_`i'^`=fp_2-1' )
		}
		if (fp_1 == fp_2) {
			gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*( fp_2*(price_`i'^`=fp_2-1')*log(price_`i') + (price_`i'^`=fp_2')*(1/price_`i') ) )
		}
	}
	qui replace rate_`i' = lower_bound if rate_`i' < lower_bound
	qui replace rate_`i' = rate_00 if rate_00 < lower_bound 
	gen emissions_`i' = rate_`i'*gas_0
	gen reduction_`i' = emissions_`=(`i'-1)' - emissions_`i'
	gen value_`i' = reduction_`i'*price_00 
	gen cost_`i' = ((price_`i' + price_`=(`i'-1)')/2)*(emissions_`=(`i'-1)' - emissions_`i')
	egen sum_cost_`i' = sum(cost_`i')
	egen sum_reduction_`i' = sum(reduction_`i')
	egen sum_value_`i' = sum(value_`i')
	gen ag_cost_`i' = ag_cost_`=`i'-1' + sum_cost_`i'
	gen ag_reduction_`i' = ag_reduction_`=`i'-1' + sum_reduction_`i'
	gen ag_value_`i' = ag_value_`=`i'-1' + sum_value_`i'
	gen mc_`i' = sum_cost_`i'/sum_reduction_`i'
	drop emissions_`=(`i'-1)' reduction_`i' value_`i' cost_`i' sum_* rate_`=(`i'-1)' price_`=(`i'-1)' //Dropping unnecessary variables to increase speed
}

do "$DATA/3_clean_sim_results.do"

preserve
keep mc_mcf mc_tco2e reduction_percent ag_reduction_tco2e
compress*
save "$DATA/Temp/sim_robust_`name'_graph.dta", replace
restore

// do "$DATA/3_sim_graph_macc.do" // Optional graph

gen check = "`name'"
keep if tax_tco2 == 5 | tax_tco2 == 20 | tax_tco2 == 44.4
keep check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
order check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
compress *
save "$DATA/Temp/sim_robust_`name'_table.dta", replace

local FE i.region_year i.facility_id // Reset

*/







// 7 // Basin-Year FE

local name basin_year_FE

use "$DATA/ghgrp_di_intermediate", clear

local FE i.basin_year i.facility_id

scalar dimension = 1
fp <p_spot>, replace dimension(`=dimension') powers(): reg `dependent' <p_spot> `controls' `FE'  `conditions', vce(cluster facility)
keep `conditions'

do "$DATA/3_sim_setup.do"

forval i = 1/`iterations' {
	gen price_`i' = price_00 + `ch4_cont'*`step_size'*`i'
	
	// First derivative
	if dimension == 1 & fp_1 != 0 {
		gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' )
	}
	if dimension == 2 & fp_1 != 0 & fp_2 != 0 {
		if (fp_1 != fp_2) {
			gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*fp_2*price_`i'^`=fp_2-1' )
		}
		if (fp_1 == fp_2) {
			gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*( fp_2*(price_`i'^`=fp_2-1')*log(price_`i') + (price_`i'^`=fp_2')*(1/price_`i') ) )
		}
	}
	qui replace rate_`i' = lower_bound if rate_`i' < lower_bound
	qui replace rate_`i' = rate_00 if rate_00 < lower_bound 
	gen emissions_`i' = rate_`i'*gas_0
	gen reduction_`i' = emissions_`=(`i'-1)' - emissions_`i'
	gen value_`i' = reduction_`i'*price_00 
	gen cost_`i' = ((price_`i' + price_`=(`i'-1)')/2)*(emissions_`=(`i'-1)' - emissions_`i')
	egen sum_cost_`i' = sum(cost_`i')
	egen sum_reduction_`i' = sum(reduction_`i')
	egen sum_value_`i' = sum(value_`i')
	gen ag_cost_`i' = ag_cost_`=`i'-1' + sum_cost_`i'
	gen ag_reduction_`i' = ag_reduction_`=`i'-1' + sum_reduction_`i'
	gen ag_value_`i' = ag_value_`=`i'-1' + sum_value_`i'
	gen mc_`i' = sum_cost_`i'/sum_reduction_`i'
	drop emissions_`=(`i'-1)' reduction_`i' value_`i' cost_`i' sum_* rate_`=(`i'-1)' price_`=(`i'-1)' //Dropping unnecessary variables to increase speed
}

do "$DATA/3_clean_sim_results.do"

preserve
keep mc_mcf mc_tco2e reduction_percent ag_reduction_tco2e
compress*
save "$DATA/Temp/sim_robust_`name'_graph.dta", replace
restore

// do "$DATA/3_sim_graph_macc.do" // Optional graph

gen check = "`name'"
keep if tax_tco2 == 5 | tax_tco2 == 20 | tax_tco2 == 44.4
keep check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
order check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
compress *
save "$DATA/Temp/sim_robust_`name'_table.dta", replace

local FE i.region_year i.facility_id // Reset

*/







// 8 // Excluding Mountain Region

local name no_mountain

use "$DATA/ghgrp_di_intermediate", clear

drop if region == 4

preserve
bysort facility_id: egen min_rate = min(rate)
bysort facility_id: egen max_rate = max(rate)
collapse min_rate max_rate, by(facility_id)
sort min_rate
gen trim_5 = 1 if _n <= round(_N/20)
sort max_rate
replace trim_5 = 1 if _n > _N - round(_N/20)
drop min max
sort facility_id
save "$DATA/Temp/trim_facilities", replace
restore
merge m:1 facility_id using "$DATA/Temp/trim_facilities"
drop _merge

scalar dimension = 2
fp <p_spot>, replace dimension(`=dimension') powers(): reg `dependent' <p_spot> `controls' `FE' `conditions', vce(cluster facility)
keep `conditions'

do "$DATA/3_sim_setup.do"

forval i = 1/`iterations' {
	gen price_`i' = price_00 + `ch4_cont'*`step_size'*`i'
	
	// First derivative
	if dimension == 1 & fp_1 != 0 {
		gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' )
	}
	if dimension == 2 & fp_1 != 0 & fp_2 != 0 {
		if (fp_1 != fp_2) {
			gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*fp_2*price_`i'^`=fp_2-1' )
		}
		if (fp_1 == fp_2) {
			gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*( fp_2*(price_`i'^`=fp_2-1')*log(price_`i') + (price_`i'^`=fp_2')*(1/price_`i') ) )
		}
	}
	qui replace rate_`i' = lower_bound if rate_`i' < lower_bound
	qui replace rate_`i' = rate_00 if rate_00 < lower_bound 
	gen emissions_`i' = rate_`i'*gas_0
	gen reduction_`i' = emissions_`=(`i'-1)' - emissions_`i'
	gen value_`i' = reduction_`i'*price_00 
	gen cost_`i' = ((price_`i' + price_`=(`i'-1)')/2)*(emissions_`=(`i'-1)' - emissions_`i')
	egen sum_cost_`i' = sum(cost_`i')
	egen sum_reduction_`i' = sum(reduction_`i')
	egen sum_value_`i' = sum(value_`i')
	gen ag_cost_`i' = ag_cost_`=`i'-1' + sum_cost_`i'
	gen ag_reduction_`i' = ag_reduction_`=`i'-1' + sum_reduction_`i'
	gen ag_value_`i' = ag_value_`=`i'-1' + sum_value_`i'
	gen mc_`i' = sum_cost_`i'/sum_reduction_`i'
	drop emissions_`=(`i'-1)' reduction_`i' value_`i' cost_`i' sum_* rate_`=(`i'-1)' price_`=(`i'-1)' //Dropping unnecessary variables to increase speed
}

do "$DATA/3_clean_sim_results.do"

preserve
keep mc_mcf mc_tco2e reduction_percent ag_reduction_tco2e
compress*
save "$DATA/Temp/sim_robust_`name'_graph.dta", replace
restore

// do "$DATA/3_sim_graph_macc.do" // Optional graph

gen check = "`name'"
keep if tax_tco2 == 5 | tax_tco2 == 20 | tax_tco2 == 44.4
keep check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
order check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
compress *
save "$DATA/Temp/sim_robust_`name'_table.dta", replace

*/






// 9 // Excluding 2016

local name no_2016

use "$DATA/ghgrp_di_intermediate", clear

drop if year == 2016

preserve
bysort facility_id: egen min_rate = min(rate)
bysort facility_id: egen max_rate = max(rate)
collapse min_rate max_rate, by(facility_id)
sort min_rate
gen trim_5 = 1 if _n <= round(_N/20)
sort max_rate
replace trim_5 = 1 if _n > _N - round(_N/20)
drop min max
sort facility_id
save "$DATA/Temp/trim_facilities", replace
restore
merge m:1 facility_id using "$DATA/Temp/trim_facilities"
drop _merge

scalar dimension = 2
fp <p_spot>, replace dimension(`=dimension') powers(): reg `dependent' <p_spot> `controls' `FE'  `conditions', vce(cluster facility)
keep `conditions'

do "$DATA/3_sim_setup.do"

forval i = 1/`iterations' {
	gen price_`i' = price_00 + `ch4_cont'*`step_size'*`i'
	
	// First derivative
	if dimension == 1 & fp_1 != 0 {
		gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' )
	}
	if dimension == 2 & fp_1 != 0 & fp_2 != 0 {
		if (fp_1 != fp_2) {
			gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*fp_2*price_`i'^`=fp_2-1' )
		}
		if (fp_1 == fp_2) {
			gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*( fp_2*(price_`i'^`=fp_2-1')*log(price_`i') + (price_`i'^`=fp_2')*(1/price_`i') ) )
		}
	}
	qui replace rate_`i' = lower_bound if rate_`i' < lower_bound
	qui replace rate_`i' = rate_00 if rate_00 < lower_bound 
	gen emissions_`i' = rate_`i'*gas_0
	gen reduction_`i' = emissions_`=(`i'-1)' - emissions_`i'
	gen value_`i' = reduction_`i'*price_00 
	gen cost_`i' = ((price_`i' + price_`=(`i'-1)')/2)*(emissions_`=(`i'-1)' - emissions_`i')
	egen sum_cost_`i' = sum(cost_`i')
	egen sum_reduction_`i' = sum(reduction_`i')
	egen sum_value_`i' = sum(value_`i')
	gen ag_cost_`i' = ag_cost_`=`i'-1' + sum_cost_`i'
	gen ag_reduction_`i' = ag_reduction_`=`i'-1' + sum_reduction_`i'
	gen ag_value_`i' = ag_value_`=`i'-1' + sum_value_`i'
	gen mc_`i' = sum_cost_`i'/sum_reduction_`i'
	drop emissions_`=(`i'-1)' reduction_`i' value_`i' cost_`i' sum_* rate_`=(`i'-1)' price_`=(`i'-1)' //Dropping unnecessary variables to increase speed
}

do "$DATA/3_clean_sim_results.do"

preserve
keep mc_mcf mc_tco2e reduction_percent ag_reduction_tco2e
compress*
save "$DATA/Temp/sim_robust_`name'_graph.dta", replace
restore

// do "$DATA/3_sim_graph_macc.do" // Optional graph

gen check = "`name'"
keep if tax_tco2 == 5 | tax_tco2 == 20 | tax_tco2 == 44.4
keep check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
order check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
compress *
save "$DATA/Temp/sim_robust_`name'_table.dta", replace

*/







// 10 // Lower bound at bottom tenth of facilities

local name lower_bound

use "$DATA/ghgrp_di_intermediate", clear

scalar dimension = 2
fp <p_spot>, replace dimension(`=dimension') powers(-2 -2): reg `dependent' <p_spot> `controls' `FE'  `conditions', vce(cluster facility)
keep `conditions'

sort rate
local lower = `=rate[107]'
di `lower'

do "$DATA/3_sim_setup.do"

replace lower_bound = `lower'

forval i = 1/`iterations' {
	gen price_`i' = price_00 + `ch4_cont'*`step_size'*`i'
	
	// First derivative
	if dimension == 1 & fp_1 != 0 {
		gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' )
	}
	if dimension == 2 & fp_1 != 0 & fp_2 != 0 {
		if (fp_1 != fp_2) {
			gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*fp_2*price_`i'^`=fp_2-1' )
		}
		if (fp_1 == fp_2) {
			gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*( fp_2*(price_`i'^`=fp_2-1')*log(price_`i') + (price_`i'^`=fp_2')*(1/price_`i') ) )
		}
	}
	qui replace rate_`i' = lower_bound if rate_`i' < lower_bound
	qui replace rate_`i' = rate_00 if rate_00 < lower_bound 
	gen emissions_`i' = rate_`i'*gas_0
	gen reduction_`i' = emissions_`=(`i'-1)' - emissions_`i'
	gen value_`i' = reduction_`i'*price_00 
	gen cost_`i' = ((price_`i' + price_`=(`i'-1)')/2)*(emissions_`=(`i'-1)' - emissions_`i')
	egen sum_cost_`i' = sum(cost_`i')
	egen sum_reduction_`i' = sum(reduction_`i')
	egen sum_value_`i' = sum(value_`i')
	gen ag_cost_`i' = ag_cost_`=`i'-1' + sum_cost_`i'
	gen ag_reduction_`i' = ag_reduction_`=`i'-1' + sum_reduction_`i'
	gen ag_value_`i' = ag_value_`=`i'-1' + sum_value_`i'
	gen mc_`i' = sum_cost_`i'/sum_reduction_`i'
	drop emissions_`=(`i'-1)' reduction_`i' value_`i' cost_`i' sum_* rate_`=(`i'-1)' price_`=(`i'-1)' //Dropping unnecessary variables to increase speed
}

do "$DATA/3_clean_sim_results.do"

preserve
keep mc_mcf mc_tco2e reduction_percent ag_reduction_tco2e
compress*
save "$DATA/Temp/sim_robust_`name'_graph.dta", replace
restore

// do "$DATA/3_sim_graph_macc.do" // Optional graph

gen check = "`name'"
keep if tax_tco2 == 5 | tax_tco2 == 20 | tax_tco2 == 44.4
keep check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
order check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
compress *
save "$DATA/Temp/sim_robust_`name'_table.dta", replace

*/






// 11 // Start facilities at 2016 values

local name start_2016

use "$DATA/ghgrp_di_intermediate", clear

scalar dimension = 2
fp <p_spot>, replace dimension(`=dimension') powers(-2 -2): reg `dependent' <p_spot> `controls' `FE'  `conditions', vce(cluster facility)
keep `conditions'

// Recovering parameters
mat coef = e(b)
mat fp_pwr =  e(fp_fp)
scalar beta_1 = coef[1,1]
scalar constant = coef[1,`=colsof(coef)']
scalar fp_1 = fp_pwr[1,1]
if dimension == 2 {
	scalar fp_2 = fp_pwr[1,2]
	scalar beta_2 = coef[1,2]
}

bysort facility_id: egen max_year = max(year) // Comment out collapse above and bring these two lines in for starting in 2016 robustness check
keep if year == max_year
keep facility_id rate p_spot gas

egen total_gas = sum(gas) // Total gas in sample
replace gas = gas*(32592000000/total_gas) // Scaling up gas production based on EIA 2016 estimate (from EIA March 2020 Monthly Energy Report)
gen emissions = rate*gas
sum emissions
scalar total_emissions = r(sum)
scalar total_emissions_tco2e = total_emissions*34/53.68

egen lower_bound = min(rate)
order lower_bound, after(facility_id)
rename rate rate_0
rename gas gas_0
rename p_spot price_0
gen emissions_0 = rate_0*gas_0
gen ag_cost_0 = 0
gen ag_reduction_0 = 0
gen ag_value_0 = 0
gen rate_00 = rate_0 // Reference so others can be deleted
gen price_00 = price_0 // Reference so others can be deleted

forval i = 1/`iterations' {
	gen price_`i' = price_00 + `ch4_cont'*`step_size'*`i'
	
	// First derivative
	if dimension == 1 & fp_1 != 0 {
		gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' )
	}
	if dimension == 2 & fp_1 != 0 & fp_2 != 0 {
		if (fp_1 != fp_2) {
			gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*fp_2*price_`i'^`=fp_2-1' )
		}
		if (fp_1 == fp_2) {
			gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*( fp_2*(price_`i'^`=fp_2-1')*log(price_`i') + (price_`i'^`=fp_2')*(1/price_`i') ) )
		}
	}	
	qui replace rate_`i' = lower_bound if rate_`i' < lower_bound
	qui replace rate_`i' = rate_00 if rate_00 < lower_bound 
	gen emissions_`i' = rate_`i'*gas_0
	gen reduction_`i' = emissions_`=(`i'-1)' - emissions_`i'
	gen value_`i' = reduction_`i'*price_00 
	gen cost_`i' = ((price_`i' + price_`=(`i'-1)')/2)*(emissions_`=(`i'-1)' - emissions_`i')
	egen sum_cost_`i' = sum(cost_`i')
	egen sum_reduction_`i' = sum(reduction_`i')
	egen sum_value_`i' = sum(value_`i')
	gen ag_cost_`i' = ag_cost_`=`i'-1' + sum_cost_`i'
	gen ag_reduction_`i' = ag_reduction_`=`i'-1' + sum_reduction_`i'
	gen ag_value_`i' = ag_value_`=`i'-1' + sum_value_`i'
	gen mc_`i' = sum_cost_`i'/sum_reduction_`i'
	drop emissions_`=(`i'-1)' reduction_`i' value_`i' cost_`i' sum_* rate_`=(`i'-1)' price_`=(`i'-1)' //Dropping unnecessary variables to increase speed
}


keep ag_* mc_*
keep if _n == 1

gen dummy = 1
reshape long ag_cost_ ag_reduction_ ag_value_ mc_, i(dummy) j(tax)
drop dummy

rename tax tax_mcf
rename mc_ mc_mcf
rename ag_reduction_ ag_reduction_mcf
rename ag_cost_ ag_cost
rename ag_value_ ag_value
replace mc_mcf = mc_mcf + (3.123 - `=mc_mcf[2]') // Adjust so facilities start at average prices for consistency with other simulations

gen marginal_value = 3.123

replace mc_mcf = marginal_value if tax_mcf == 0

replace tax_mcf = tax_mcf*$step_size
recast double tax_mcf
replace tax_mcf = round(tax_mcf,.01)
gen tax_tco2 = (tax*53.68)/34
gen ag_reduction_tco2e = ag_reduction_mcf*34/53.68
gen mc_tco2e = mc_mcf*53.68/34
di total_emissions
gen reduction_percent = (ag_reduction_mcf/total_emissions)*100

gen net_cost = (ag_cost - ag_value)/28400000000

replace ag_reduction_mcf = 0 if missing(ag_reduction_mcf)
replace ag_reduction_tco2e = 0 if missing(ag_reduction_tco2e)

format mc* tax* %6.2fc
format ag* %12.0fc
format reduction_percent %6.1fc

order tax_mcf tax_tco2 mc_mcf mc_tco2e ag_cost ag_reduction_mcf ag_reduction_tco2e reduction_percent

// Interpolating $5, $20, and $42 carbon prices
set obs `=_N+1'
replace tax_tco2 = 5 if missing(tax_tco2)
set obs `=_N+1'
replace tax_tco2 = 20 if missing(tax_tco2)
set obs `=_N+1'
replace tax_tco2 = 44.4 if missing(tax_tco2)
sort tax_tco2
gen percent_next = (tax_tco2 - tax_tco2[_n-1])/(tax_tco2[_n+1] - tax_tco2[_n-1]) if missing(mc_mcf)
gen percent_prev = (tax_tco2[_n+1] - tax_tco2)/(tax_tco2[_n+1] - tax_tco2[_n-1]) if missing(mc_mcf)
replace tax_mcf = tax_mcf[_n-1]*percent_prev + tax_mcf[_n+1]*percent_next if missing(tax_mcf)
replace tax_tco2 = tax_tco2[_n-1]*percent_prev + tax_tco2[_n+1]*percent_next if missing(tax_tco2)
foreach var in mc_mcf mc_tco2e ag_cost ag_reduction_mcf ag_reduction_tco2e reduction_percent ag_value net_cost {
	replace `var' = `var'[_n-1]*percent_prev + `var'[_n+1]*percent_next if missing(`var')
}

gen avg_cost_tco2e = (ag_cost - ag_value)/ag_reduction_tco2e
format avg_cost_tco2e %6.2fc
format net_cost %6.4fc

recast double tax_tco2
replace tax_tco2 = 44.4 if tax_tco2 > 44.39 & tax_tco2 < 44.41 // Became 44.400002 due to some variable formatting bug

gen double ag_reduction_tco2e_round = round(ag_reduction_tco2e,1000)
format ag_reduction_tco2e_round %15.0fc
drop ag_reduction_tco2e
rename ag_reduction_tco2e_round ag_reduction_tco2e

gen total_net_cost = (ag_cost - ag_value)/1000000
gen climate_benefit = (ag_reduction_tco2e*44.4)/1000000
gen welfare = climate_benefit - total_net_cost

format climate_benefit welfare total_net_cost %10.0fc

replace ag_cost = ag_cost/1000000
replace ag_value = ag_value/1000000

preserve
keep mc_mcf mc_tco2e reduction_percent ag_reduction_tco2e
compress*
save "$DATA/Temp/sim_robust_`name'_graph.dta", replace
restore

// do "$DATA/3_sim_graph_macc.do" // Optional graph (note y-axis issue with this one; doesn't matter for results below)

gen check = "`name'"
keep if tax_tco2 == 5 | tax_tco2 == 20 | tax_tco2 == 44.4
keep check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
order check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
compress *
save "$DATA/Temp/sim_robust_`name'_table.dta", replace

*/






// 12 // Modelling facility size - this one requires a little more detail

local name size

use "$DATA/ghgrp_di_intermediate", clear

keep if trim_5 != 1

gen p_spot_gas = p_spot*mean_gas
gen p_gas_fp1 = p_spot_gas^(-.5)
gen p_gas_fp2 = p_spot_gas^(-.5)*log(p_spot_gas)

// Second-order FP Regression
eststo: fp <p_spot>, replace dimension(2) powers(): reg `dependent' <p_spot> p_gas_fp1 p_gas_fp2 `controls' `FE'  `conditions', vce(cluster facility) // FP coefficients are -2 -1

// Recovering parameters
mat coef = e(b)
mat fp_pwr =  e(fp_fp)
scalar beta_1 = coef[1,1]
scalar beta_2 = coef[1,2]
scalar beta_3 = coef[1,3]
scalar beta_4 = coef[1,4]
scalar constant = coef[1,`=colsof(coef)']
scalar fp_1 = fp_pwr[1,1]
scalar fp_2 = fp_pwr[1,2]

collapse (mean) rate gas p_spot, by(facility_id)

gen mean_gas = gas
gen p_spot_gas = p_spot*mean_gas

egen total_gas = sum(gas) // Total gas in sample
replace gas = gas*(32592000000/total_gas) // Scaling up gas production based on EIA 2016 estimate (from EIA March 2020 Monthly Energy Report)
gen emissions = rate*gas
sum emissions
scalar total_emissions = r(sum)
scalar total_emissions_tco2e = total_emissions*34/53.68

egen lower_bound = min(rate)
order lower_bound, after(facility_id)
rename rate rate_0
rename gas gas_0
rename p_spot price_0
rename p_spot_gas price_gas_0 // Added
gen emissions_0 = rate_0*gas_0
gen ag_cost_0 = 0
gen ag_reduction_0 = 0
gen ag_value_0 = 0
gen rate_00 = rate_0 // Reference so others can be deleted
gen price_00 = price_0 // Reference so others can be deleted
gen price_gas_00 = price_gas_0 // Reference so others can be deleted

forval i = 1/`iterations' {
	gen price_`i' = price_00 + `ch4_cont'*`step_size'*`i'
	gen price_gas_`i' = price_gas_00 + mean_gas*`ch4_cont'*`step_size'*`i'
	gen rate_`i'= rate_`=(`i'-1)' + `ch4_cont'*`step_size'*(fp_1*beta_1*price_`i'^`=fp_1-1' + beta_2*(fp_2*price_`i'^`=fp_2-1') + beta_3*(-.5*price_gas_`i'^(-1.5)) + beta_4*( (-.5*price_gas_`i'^(-1.5))*log(price_gas_`i') + (price_gas_`i'^(-.5))*(1/price_gas_`i') ) ) // First derivative
	qui replace rate_`i' = lower_bound if rate_`i' < lower_bound
	qui replace rate_`i' = rate_00 if rate_00 < lower_bound // Reference start point
	
	gen emissions_`i' = rate_`i'*gas_0
	gen reduction_`i' = emissions_`=(`i'-1)' - emissions_`i'
	gen value_`i' = reduction_`i'*price_00 // Value of conserved gas
	gen cost_`i' = ((price_`i' + price_`=(`i'-1)')/2)*(emissions_`=(`i'-1)' - emissions_`i')
	
	egen sum_cost_`i' = sum(cost_`i')
	egen sum_reduction_`i' = sum(reduction_`i')
	egen sum_value_`i' = sum(value_`i')
	gen ag_cost_`i' = ag_cost_`=`i'-1' + sum_cost_`i'
	gen ag_reduction_`i' = ag_reduction_`=`i'-1' + sum_reduction_`i'
	gen ag_value_`i' = ag_value_`=`i'-1' + sum_value_`i'
	gen mc_`i' = sum_cost_`i'/sum_reduction_`i'
	
	drop emissions_`=(`i'-1)' reduction_`i' value_`i' cost_`i' sum_* rate_`=(`i'-1)' price_`=(`i'-1)' //Dropping unnecessary variables to increase speed
}

do "$DATA/3_clean_sim_results.do"

preserve
keep mc_mcf mc_tco2e reduction_percent ag_reduction_tco2e
compress*
save "$DATA/Temp/sim_robust_`name'_graph.dta", replace
restore

do "$DATA/3_sim_graph_macc.do" // Optional graph
graph export "$DATA/Graphs/p8_macc_size.pdf", replace

gen check = "`name'"
keep if tax_tco2 == 5 | tax_tco2 == 20 | tax_tco2 == 44.4
keep check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
order check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
compress *
save "$DATA/Temp/sim_robust_`name'_table.dta", replace

*/




// 13 // Negative Binomial Model

local name NB

use "$DATA/ghgrp_di_intermediate", clear

scalar dimension = 2
eststo: fp <`independent'>, dimension(2) powers(-2 -1) replace: nbreg ch4_normal <`independent'> `controls' `FE' `conditions', exposure(gas) vce(cluster facility) // 2nd-degree FP
keep `conditions'

// Recovering parameters
mat coef = e(b)
mat fp_pwr =  e(fp_fp)
scalar beta_1 = coef[1,1]
scalar constant = coef[1,`=colsof(coef)']
scalar fp_1 = fp_pwr[1,1]
if dimension == 2 {
	scalar fp_2 = fp_pwr[1,2]
	scalar beta_2 = coef[1,2]
}

collapse (mean) rate gas p_spot, by(facility_id)

egen total_gas = sum(gas) // Total gas in sample
replace gas = gas*(32592000000/total_gas) // Scaling up gas production based on EIA 2016 estimate (from EIA March 2020 Monthly Energy Report)
gen emissions = rate*gas
sum emissions
scalar total_emissions = r(sum)
scalar total_emissions_tco2e = total_emissions*34/53.68

gen log_rate = log(rate)
egen lower_bound = min(log_rate)
order lower_bound, after(facility_id)

rename log_rate log_rate_0
rename rate rate_0
rename gas gas_0
rename p_spot price_0
gen emissions_0 = rate_0*gas_0
gen ag_cost_0 = 0
gen ag_reduction_0 = 0
gen ag_value_0 = 0
gen log_rate_00 = log_rate_0 // Reference so others can be deleted
gen rate_00 = rate_0 // Reference so others can be deleted
gen price_00 = price_0 // Reference so others can be deleted



forval i = 1/`iterations' {
	gen price_`i' = price_00 + `ch4_cont'*`step_size'*`i'
	
	// First derivative
	gen log_rate_`i'= log_rate_`=(`i'-1)' + `ch4_cont'*`step_size'*( beta_1*fp_1*price_`i'^`=fp_1-1' + beta_2*fp_2*price_`i'^`=fp_2-1' ) // Specified for NB, FP coefficients not equal
	
	qui replace log_rate_`i' = lower_bound if log_rate_`i' < lower_bound
	qui replace log_rate_`i' = log_rate_00 if log_rate_00 < lower_bound
	gen rate_`i' = exp(log_rate_`i')
	
	gen emissions_`i' = rate_`i'*gas_0
	gen reduction_`i' = emissions_`=(`i'-1)' - emissions_`i'
	gen value_`i' = reduction_`i'*price_00 
	gen cost_`i' = ((price_`i' + price_`=(`i'-1)')/2)*(emissions_`=(`i'-1)' - emissions_`i')
	egen sum_cost_`i' = sum(cost_`i')
	egen sum_reduction_`i' = sum(reduction_`i')
	egen sum_value_`i' = sum(value_`i')
	gen ag_cost_`i' = ag_cost_`=`i'-1' + sum_cost_`i'
	gen ag_reduction_`i' = ag_reduction_`=`i'-1' + sum_reduction_`i'
	gen ag_value_`i' = ag_value_`=`i'-1' + sum_value_`i'
	gen mc_`i' = sum_cost_`i'/sum_reduction_`i'
	drop emissions_`=(`i'-1)' reduction_`i' value_`i' cost_`i' sum_* rate_`=(`i'-1)' price_`=(`i'-1)' //Dropping unnecessary variables to increase speed
}

do "$DATA/3_clean_sim_results.do"

preserve
keep mc_mcf mc_tco2e reduction_percent ag_reduction_tco2e
compress*
save "$DATA/Temp/sim_robust_`name'_graph.dta", replace
restore

// do "$DATA/3_sim_graph_macc.do" // Optional graph

gen check = "`name'"
keep if tax_tco2 == 5 | tax_tco2 == 20 | tax_tco2 == 44.4
keep check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
order check tax_mcf tax_tco2 reduction_percent ag_reduction_tco2e ag_cost ag_value total_net_cost net_cost climate_benefit welfare
compress *
save "$DATA/Temp/sim_robust_`name'_table.dta", replace

*/







// 14 // Combine

// 14.1 // Graph

use "$DATA/Temp/sim_robust_main_graph.dta", clear
gen check = "main"
foreach name in FP1 weighted trim1 state_regs comp_lastyr year_fe basin_year_fe no_mountain no_2016 size lower_bound start_2016 NB {
	append using "$DATA/Temp/sim_robust_`name'_graph.dta"
	replace check = "`name'" if missing(check)
}
drop ag_reduction_tco2e
drop mc_tco2e
gen marginal_value = 3.123

// Graph

// Value of conserved gas for area plot
set obs `=_N+1'
replace marginal_value = marginal_value[_n-1] if missing(marginal_value)
replace reduction_percent = 100 if missing(reduction_percent)

// Indicating out-of-sample predictions: Only 79% that is methane is taxed, so the highest observe price (5.75) is equivalent to a methane tax of 5.75/.79 = 7.278
set obs `=_N+1'
replace reduction_percent = 98.75 if _n == _N
gen reduction_percent_end = 98.75 if _n == _N
replace mc_mcf = 7.278 if _n == _N
gen mc_mcf_end = 29.75 if _n == _N

set obs `=_N+2' // This makes a little dash
replace mc_mcf = 7.278 if _n >= _N - 1
replace reduction_percent = 98 if _n == _N - 1
replace reduction_percent = 99.5 if _n == _N

twoway area marginal_value reduction_percent,  color(gs14) lwidth(vvthin) lcolor(gs14) || ///
	line mc_mcf reduction_percent if check == "FP1", 			lcolor(black) 		lpattern(longdash) 	lwidth(medthick) || ///
	line mc_mcf reduction_percent if check == "year_fe", 		lcolor(black) 		lpattern(#-#)		lwidth(medthick) || ///
	line mc_mcf reduction_percent if check == "basin_year_fe", 	lcolor(black) 		lpattern(.#) 		lwidth(medthick) 	 || ///
	line mc_mcf reduction_percent if check == "NB", 			lcolor(black) 		lpattern() 			lwidth(thin) || ///
	line mc_mcf reduction_percent if check == "state_regs", 	lcolor(ebblue) 							lwidth(medthick) || ///
	line mc_mcf reduction_percent if check == "comp_lastyr", 	lcolor(ebblue) 		lpattern(longdash) 	lwidth(medthick) || ///
	line mc_mcf reduction_percent if check == "no_mountain", 	lcolor(ebblue) 		lpattern(#-#) 		lwidth(medthick) || ///
	line mc_mcf reduction_percent if check == "no_2016", 		lcolor(ebblue) 		lpattern(.#) 		lwidth(medthick) 	 || ///
	line mc_mcf reduction_percent if check == "trim1", 			lcolor(ebblue) 	 						lwidth(thin) || ///
	line mc_mcf reduction_percent if check == "size", 			lcolor(olive_teal) 	 					lwidth(medthick) || ///
	line mc_mcf reduction_percent if check == "weighted", 		lcolor(olive_teal) 	lpattern(longdash)	lwidth(medthick) || ///
	line mc_mcf reduction_percent if check == "lower_bound", 	lcolor(olive_teal) 	lpattern(#-#) 		lwidth(medthick) || ///
	line mc_mcf reduction_percent if check == "start_2016", 	lcolor(olive_teal) 	lpattern(.#) 		lwidth(medthick) 	 || ///
	line mc_mcf reduction_percent if check == "main", 			lcolor(black) 							lwidth(thick) 	 || ///
	line mc_mcf reduction_percent if _n >= _N - 1, lcolor(gs12) lwidth(thin) || ///
	pcarrow mc_mcf reduction_percent mc_mcf_end reduction_percent_end if _n == _N - 2, lcolor(gs12) mcolor(gs12) mlwidth(thin) msize(1.75) barbsize(.35) lwidth(thin) ///
	yline(5.68, lcolor(gs12) lwidth(vthin) ) /// Obtain these by examining mc_mcf in the dataset at i.e. tax_tco2 = 5
	yline(12.67, lcolor(gs12) lwidth(vthin) ) ///
	yline(25.41, lcolor(gs12) lwidth(vthin) ) ///
	text(5.68 .5 "$5 Carbon Price" , yaxis(1) xaxis(1) placement(ne) orient(hor) margin(.5 .5 .5 .5) color(gs8) size(medsmall) ) ///
	text(12.67 .5 "$20 Carbon Price", yaxis(1) xaxis(1) placement(ne) orient(hor) margin(.5 .5 .5 .5) color(gs8) size(medsmall) ) ///
	text(25.41 .5 "Social Cost of Methane", yaxis(1) xaxis(1) placement(ne) orient(hor) margin(.5 .5 .5 .5) color(gs8) size(medsmall) ) ///
	text(.5 50 "Value of Captured Gas", yaxis(1) xaxis(1) placement(n) orient(hor) margin(.5 .5 .5 .5) color(gs8) size(small) ) ///
	text(`=(29.75+7.278)/2' 100 "Out-of-Sample Predictions" , yaxis(1) xaxis(1) placement(e) orient(vert) margin(.5 .5 .5 .5) color(gs12) size(*.6)) ///
	graphregion(color(white) margin(3 4 2 0)) xsize(7) ysize(4) ///
	plotregion(margin(0 0 0 0)) ///
	xlabel(0 "0%" 20 "20%" 40 "40%" 60 "60%" 80 "80%" 100 "100%", axis(1) labsize(medsmall) format(%12.0fc) nogrid ) ///
	xscale(range(0 100)) ///
	xtitle("Total Abatement"  , axis(1) size() margin(0 0 22 28)) ///
	title(" "  , margin(0 0 0 5)) ///
	ylabel(0(5)31, axis(1) labsize(medsmall)format(%6.0fc) nogrid) ///
	yscale(range(0 31)) ///
	ytitle("Marginal Cost ($/Mcf)", axis(1) size() margin(0 2 0 0)) ///
	legend(order(15 "Main Specification" 2 "First-Order FP" 3 "Year FE" 4 "Basin-Year FE" 5 "Negative" "Binomial Model" 6 "Trimming Based" "On State Regs." 7 "Controlling for Prev." "Year Completions" 8 "Removing Mountain" "Region" ///
		9 "Removing 2016" 10 "Trimming Emission" "Rates at 1%" 11 "Heterogeneous" "Facility Size" 12 "Weighting by Avg" "Gas Production" 13 "Lower Bound at" "Bottom Decile" 14 "Starting at 2016" "Prices & Rates" ) ///
		ring(2) size(*.8) linegap(.25) rowgap(3.5) region(lwidth(thin) lcolor(gs14)  margin(2 1 1 1)) symxsize(7.25) margin(1 1 1 1) bmargin(12 0 -35 -5) pos(3) col(1))
graph export "$DATA/Results/p8_sim_robust.png", width(2800) height(1600) replace

*/

// 13.2 // Table
use "$DATA/Temp/sim_robust_main_table.dta", clear
foreach name in FP1 weighted trim1 state_regs comp_lastyr year_fe basin_year_fe no_mountain no_2016 size lower_bound start_2016 NB {
	append using "$DATA/Temp/sim_robust_`name'_table.dta"
}
drop if tax_tco2 == 20
sort tax_mcf reduction_percent

export excel "$DATA/Results/sim_robust.xls", firstrow(var) replace
